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ABSTRACT 

In data from three clear nights of a WHT/UES run in 2000 Oct /Nov, and using im- 
proved Doppler tomographic signal-analysis techniques, we have carried out a deep 
search for starlight reflected from the innermost of v (upsilon) And's three plan- 
ets. We place upper limits on the planet's radius Rp as functions of its projected 
orbital velocity Kp « 139 sin z km s^^ for various assumptions about the wavelength- 
dependent geometric albedo spectrum p(A) of its atmosphere. For a g rey albedo p we 



find Rp^/v < 0.98 R j with 0.1% false-alarm probability (4-<t). For a Sudarsky, Bur- 



rows fc Pinto (1999)1 Class V model atmosphere, the mean albedo in our 380-676 nm 



bandpass is (p) ~ 0.42, requiring Rp < 1.51 Rj, while an (isolated) Class IV model 
with (p) ~ 0.19 requires Rp < 2.23 Rj. The star's Vrot sini ^ 10 km s^^ and estimated 
rotation period Prot lOd suggest a high orbital inchnation i ~ 70 — 80°. We also 
develop methods for assessing the false-alarm probabilities of faint candidate detec- 
tions, and for extracting information about the albedo spectrum and other planetary 
parameters from faint reflected-light signals. 
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1 INTRODUCTION 

The discovery of "hot Jupiters" - giant exoplanets in 3 
to 5-dav orbits about their parent stars ( Mayor fc Queloz 



1995 ; Marcy fc Butler 1998 ) - has overturned conventional 
theories of planetary-system formation. Theorists working 
to explain the solar system had posited that Jupiters could 
form only outside 5 AU, where planetesimals retain ice man- 
tles. Inward migration was considered possible, but deemed 



(1) the age at which they arrived in their present close or- 
bits, and (2) their atmospheric albedos, which control the 
atmospheric ter nperature-pressure g radient and hence the 
rate of cooling ( Burrows et al. 200Cl| ). 

Th e discovery of transits o f the planet orbiting HD 
209458 ( |Charbonneau et al. 200C ) provided conclusive proof 
of the existence and gas-giant nature of the "hot Jupiters" . 
Ultra-precise HST/STIS photometry of HD 20 9458 b's tran 



unlikelji since the planets would then be swallowed up hy 



sit p rofile yields a radius about 1.35 Rjup (Brown et al 



the par ent otar. Hot Jupitoro have now revived the idea of 
inward migration, but require mechanisms that halt some of 
the pla iiGto near 0.05 AU. 



2001), although the planet has a mass only 0.63 Mjup- The 
next crucial step will be to measure albedos directly. 



Theoretical studies of the structure and e^^olution of ea- 
oplanets are yielding testable predictions. Gas-giant models 
yield planet radii Rj, ~ 1.0—1.7 Rjuv for 1 < M^/Mjuy < 10 
(saumon et al. 1996; Guillot et al. 1996| ; Burrows et al 



S pectral models(3eager fc Sasselov 1998 ; Marley et al 



199£ ; Sudarsky, Burrows fc Pinto 1999[ ) include thermal 



1997). Hot Jupiters have fuUy-convective interiors, so the 



core entropy fixes the radius at a given mass ( Burrows et al 
200C). Isolated Jupiters will therefore cool and contract to 
about Jupiter's size in a few Gyr, but strongly irradiated 
giant planets find it harder to cool. Consequently the mass- 
radius-age relation for hot Jupiters is determined largely by 
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emission in the infrared and scattered starlight in the op- 
tical. The albedo is very sensitive to the depth of cloud 
formation. At the expected temperatures (Tcff ~ 1400K) 
and gravities (g > 36 m s~'^) of massive planets such as 
r Boo b, possible cloud condensates include upper cloud 
decks of MgSiOa (enstatite) and AI2O3 (corundum) at pres- 
sures ~ 5 to 10 bar, and a deeper layer of Fe grains at 
pressures a factor 2 or so higher. If silicate cloud decks form 
this deep in the atmosphere or are absent, much of the inci- 
dent starhght is absorbed in pressure-broadened absorption 
lines of neutral sodium and potassium, leaving only a weak 
Rayleigh-scattered reflection signature at blue wavelengths. 
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Figure 1. Geometric albedo spectra for the Class V (solid line), 
iso lated (dashed) and irradiated (dot-d ashed) Class IV models 
of Sudarsky, Burrows & Pinto (1999), plotted over the wave- 



length range we observed. Together with a grey model having 
p = 0.42, the Class V and isolated Class IV models were used 
to probe the wavelength dependence of candidate reflected-light 
signals. The absorption features near 0.404 and 0.589 /^m are 
due to potassium and sodium respectively. The geometric albedo 
spectrum of Jupiter (short dashes) is shown for comparison, in- 
cluding methane absorption bands at 0.54 and 0.62 /xm. The lower 
panel shows the wavelength sensitivity of our detection method, 
expressed as the deficit of detected photons Nx per unit wave- 
length due to known absorption lines in the stellar spectrum. A'^;^ 
is expressed in units of 10^^ photons per /^m, and gives the total 
photon deficit for a typical group of 4 exposures. 



Sudarsky, Burrows & Pinto (1999) predict this "Class IV" 
albedo spectrum for planets with high surface gravities. At 
the lower surface gravities <; ~ 10 m s~^ expected of lower- 
mass planets like v And b and HD 209458 b, "Class V" 
models predict silicate cloud decks forming higher in the at- 
mosphere at pressures around 0.3 bar, giving a geometric 
albedo about 60% that of Jupiter throughout most of the 
optical spectrum (Fig.^). 

Recent attempts to detect starlight reflected from 
T Boo b have produced deep upp er limits on the planet's 
albedo. Charbonneau et al. (1999) established an upper limit 

-4 



on the fully-illuminated plane t /star flux ratio eo — 10 
Collier Cameron et al. (1999) claimed a detection with an 



average eo = 7.5 x 10~° from data secured in 1998 and 
1999. The radius inferred for the planet was, however, im- 
plausibly large. Deeper obse rvations made early in 2000 
(Collier Cameron et al. 2001) did not confirm this candi- 
date detection, but produced a more stringent upper limit, 
eo < 3.5 X 10"^ between 387.4 and 586.3 nm. The appar- 



ent tidal synchronisation of r Boo with the planet's orbit 
suggests a low inclination of order 40°. At this inclination, 
the combined 1998-2000 WHT observations yielded a 99.9% 
upper limit on the geometric albedo p < 0.22, assuming a 



planet radius 1.2Rjup as predicted by Burrows et al. (2000) 



This appears to rule out the presence of a high cloud deck 
with the high reflectivity of the Class V models, which would 
give p ~ 0.4 over most of the visible spectrum. 

Here we report results from a programme of obser- 
vations designed to detect the spectroscopic signature of 
starlight reflected from the innermost of t he three giant 
planets orbiting v And ( Butler et al. 1999| ). In Section ^ 
we use the measured system parameters to determine the 
a priori probabilities for the planet's orbital velocity ampli- 
tude and t he f raction of the star's light that it intercepts. 



In Section 2.4 we discuss the expected albedo and phase 



function that the planet might display. In Sections ^ and 
^ we describe the acquisition and extraction of the echelle 
spectra. The methods used to remove the direct-starlight 
signature from the data, and to extract the planetary signal 
by combining the proflles of the thousands of photospheric 
lines recorded in each echellogram, are presented in Section 
^ A matched-filter method for measuring the strength of 
the refiected-light signal is described in Section Finally in 
Section ^ we determine upper limits on the planet's radius 
for three different models of the planet's albedo spectrum, 
and discuss the plausibility of a candidate reflected-light sig- 
nal that appears in the data at the most probable velocity 
amplitude and signal strength. 



2 SYSTEM PARAMETERS 

V And (HR 458, HD 9826) is a late-F main-sequence star 
with parameters as listed in Table 0. High-precision radial- 
yelocity studies hav e shown it to have a system of 3 planets 
(Butler et al. 199t). For the purposes of this paper we are 
concerned only with the innermost of these planets, whose 
properties (as determined directly from radial-velocity stud- 
ies or inferred using the estimated stellar parameters) are 
also summarised in Table 0. 

Starlight reflected from the orbiting planet's atmo- 
sphere leaves a detectable signature in observed spectra in 
the form of faint copies of each of the stellar absorption lines, 
Doppler shifted by the planet's orbital motion, and greatly 
diminished in brightness because the planet intercepts only 



fraction ( 1/4) (i?p /a) 



10 of the starlight and only 



part of that is reflected back into space. 

By detecting and characterizing the planetary reflected 
light signal, we in effect observe the planet/star flux ratio 
eps = fp/ f* as a function of orbital phase <j!> and wavelength 
A. The information we aim to obtain (in order of increasing 
difficulty) comprises: 

(i) Kp, the planet's projected orbital velocity. From this 
we learn the orbital inclination i and hence the planet mass 
Mp, since Mp sini is known from the star's Doppler wobble. 

(ii) eo, the maximum strength of the reflected starlight. 
From this we constrain the planet's radius since eo — 
Rp/a^fp, where p is the geometric albedo. 

(iii) p(A), the albedo spectrum, which depends on the 
composition and structure of the planetary atmosphere. 
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Table 1. Physical parameters for vAnd and its innermost planet. 



Parameter 


Value (Uncertainty) 


Star: 




Spectral type 


F8V 


Teff 


6100 (80) 


M. (Mq) 


1.3 (0.1) 


R* (Rq) 


1.63 (0.06) 


logs 


4.0 


[Fe/H] 


0.09 to 0.17 


TT (mas) 


74.25 (0.70) 


vsini (km s~^) 


9.5 (0.8) 


log Prot (days) 


1.08 (0.08) 


Age (Gyr) 


2.6 to 4.8 


Planet: 




Orbital period Port (days) 


4.61707 (0.00003) 


Transit epoch To (JD) 


2451821.457 (0.012) 


K^, (m s-i) 


74.5 (2.3) 


a (AU) 


0.059 (0.002) 


Mp sin i (Mjjjp) 


0.73 (0.04) 



Gonzalez a997) 


Fuhrmann, Pfeiffer & Bernkoof (1998): 


Gonzalez (1998) 


Gonzalez (1997) 


Fuhrmann. Pfeiffer & 


Bernkopf (1998): 


Gonzalez (1998) 


Fuhrmaun. Pfcif 


cr & Bernkopf (1998): 


Ford. Rasio k- Sills (1999' 




Fuhrmaun. PfcifFcr & Bernkoof ('1998): 


Ford. Rasio & Sills (1999' 


Gonzalez (1997) 


Fuhrmann, Pfeiffer & Bernkopf ( 1998) : 


Gonzalez (1998) 


Gonzalez (1997) 


Fuhrmann, Pfeiffer & Bernkopf (1998): 


Gonzalez (1998) 


Ferryman et al. 


'1997)1 






Fiihrmann, Pfeiffer fy. Rernkopf (1998) 






Henrv et al. (2000) BaUunas et al. (1997) Noves et al. (1984) 


Fuhrmann, Pfeiffer & Bernkopf (1998); 


Ford, Rasio & Sills (1999; 



Marcy, personal communication 
Marcv. personal com munication 



Butler et al. (19991 



Butler et al. (1999)| (revised value derived in this paper) 



Butler et al. (1999]| (revised value derived in this paper) 



(iv) g{a,X), the phase function describing the dependence 
of the amount of Ught reflected toward the observer on the 
star-planet-observer angle a. 

(v) A, the velocity width of lines in the reflected starlight. 
This depends mainly on the star's rotation in the frame of 
the orbiting planet, and to a lesser extent on the planet's 
rotation and surface winds. 

At this stage of the subject, our primary goal is to achieve 
detections of the planetary reflected light signal, thereby 
reaching step (2) in the above list - measuring Kp and eo 
and hence Rp for an assumed albedo p. With sufficient data 
we can measure epsilon{\) and hence p(A), but with present 
data the best we can do is to adopt p(A) for various plane- 
tary atmospheric models and for each one measure or place 
an upper limit on the planet radius Rp/a. To optimize the 
sensitivity for detection, we restrict observations to gibbous 
phases of the planetary orbit, thus sacrificing information 
about the phase function g{a,X). We elaborate these issues 
in more detail below. 



2.1 Radial velocity amplitude 

The radial velocity curve of the reflected light at orbital 
phase <j) is 



Vp{(j>) — KpSm2-K(j>, 



(1) 



The orbital phase at time t is given hy (j> = {t — To) /Port, 
using the values of Port and To given in Table ^. 

The mass ratio q = Mp/Mi, is determined from the 
observed amplitude Ki, of the star's reflex orbital velocity 



via; 



1 + g 



Porb 

sin i 27ra 
5.339 X 10-"* 



V74.5 m s-iy \ 1.3M0 



-1/3 



(2) 



The apparent radial velocity amplitude Kp of the planet 
about the system's centre of mass is 



27ra sini sini 
Kp = = 139 



" Poril+q (l+g)2/3 \^1.3M0 

Kepler's third law gives 

1/3 



1/3 



a = 0.0592 



1.3Mp 



AU. 



km s ^. (3) 



(4) 



The orbital inclination i is, according to the usual con- 
vention, the angle between the orbital angular momentum 
vector and the line of sight. For all but the lowest inclina- 
tions, the planet's orbital velocity amplitude is considerably 
greater than the typical widths of the star's photospheric 
absorption lines. Lines in the planet's reflected-light spec- 
trum should therefore be Doppler-shifted well clear of their 
stellar counterparts, allowing clean spectral separation for 
most of the orbit. 



2.2 Orbital inclination 

Orbital inclinations i > 80° or so are ruled out by the ab- 



sence of transits in high-precision photometry (Henry et al 
200C). The minimum inclination at which grazing transits 
occur is given by 

R* -\- Rp 



cos %rnin — 



(5) 



Low inclinations are also ruled out by the star's chromo- 
spheric activity levels. The stellar radius and Dsini suggest 
an axial rotation period of 9 days. A 12 ± 2 day period was 
estimated by Baliunas et al. (1997) based on four Ca ll H 
and K flux measurements and the period-activity-colour re- 
lat ion of Noyes et al. ( 1984). A more comprehensive study 



by Henry et al. (2000), based on an additional 212 Ca ll H 
and K flux measurements in the 1996, 1997 and 1998 ob- 
serving seasons shows the long-term average Ca ll H and K 
flux level to be almost identical to the original estimate by 
Baliunas et al. The rms scatter in the 30-day means is 4.7% 
of the long-term mean flux level. The 20% uncertainty in the 
11.7-day period estimated from the long-term mean chromo- 
spheric flux is therefore dominated by the intrinsic scatter 
in of about 0.08 dex in the Noyes et al. (1984) calibration. 
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Figure 2. The greyscale shows the prior joint probability density 
function (PDF) for projected orbital velocity Kp and the squared 
ratio (Rp/a)^ of the planet radius to the orbit radius, based on the 
measured system parameters. Darker shades in the greyscale de- 
note greater probabilities of the planet having the corresponding 
combination of (Rp/a)^ and Kp. The one-dimensional projections 
of the PDF on to the two axes are shown as histograms. They 
show that the planet is most likely to have Kp ~ 135 km s~^ and 
{Rp/a)'^ ~ 10"'*. The horizontal line shows the value of (Rp/a)^ 
for a 1 Rjup planet in the same orbit. 



Significantly faster axial rotation is needed to match the ob- 
served V sin i at orbital inclinations less than 60° , but is hard 
to reconcile with the observed chromospheric activity level. 
Henry et al. also point out that frequency analyses of the 
short-term variability in the 1996 and 1998 seasons yielded 
possible low-amplitude rotational modulation signals with 
periods of 11 and 19 days respectively, although both were 
classified as "weak". 

These rather loose constraints on i require a Monte 
Carlo simulation to express our knowledge of the system 
parameters. For this we first generate a distribution of sin i 
values 
. . Pcai I. V sin i 



(6) 



computed from randomly-chosen values of the stellar rota- 
tion period, radius and vsini and rejecting those combina- 
tions that yielded sini > 1. We assume gaussian distribu- 
tions for the log of the rotation period, the measured stellar 
mass and radius, and the stellar reflex velocity K^. 

High inclinations are then eliminated by the observed 
absence of eclipses. Low i are then given a reduced proba- 
bility consistent with the measured Vrotsini = 9.5 km s^^ 
and estimate of its rotation period Prot ~ 12 ± 2d, which 
together favour a relatively high i. The lower histogram in 
Fig. ^ shows the resulting probability distribution for Kp 
(equivalent to i). 
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Figure 3. A Monte Carlo simulation shows how the apparent 
stellar equatorial rotation speed i^e.re/i in the starlight incident on 
the planet is expected to vary with the projected orbital velocit 
amplitude Kp. The system parameters are defined as in Section 
The median rotational broadening of the reflected starlight is 8.4 
km s"-*^, somewhat less than the v sini = 9.5 km of the direct 
starlight. 



2.3 Rotational broadening 

The matched-filter analysis used to search for the reflected- 
light signal requires an a priori estimate of the rotational 
broadening of the spectral-line proflles in starlight reflected 
from the planet's atmosphere. 

The Earth-bound observer sees stellar aborption lines 
that are rotationally broadened by 



Vi, sm I = 



27ri?* sin i 



= 9.5 ± 0.8km s" 



(7) 



We assume that the planet orbits in the same direction as 
the stellar rotation and that the orbital plane and stellar 
equator are coincident. The light received by the orbiting 
planet, and subsequently reflected toward an observer in any 
direction, then has a rotational broadening 



Urefl 



27ri?* 



Prot 



1 

Porh 



(8) 



The Monte Carlo simulation gives the resulting relationship 
between the broadeniiig Vrefi of the reflected starlight and 
Kp, as shown in Fig. 0. The distribution of the predicted 
rotational broadening of the reflected starlight has a mean 
8.3 km s~ , a = 1.1 km s~ , and median 8.4 km s~^. The 
planet line widths are reduced at lower orbital inclinations, 
because the difference between the orbital and stellar rota- 
tion frequencies is reduced. 
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Figure 4. Brightness variation of model planet with orbital phase 
assuming a Venus-like phase function (solid) or a Lambert sphere 
(dashed) for orbit inclinations i = 0, 30, 60, 90° . In both cases, 
flux is measured relative to the value for illumination phase angle 
= 0. The horizontal lines show the i = 0° case, while i = 90° 
gives the greatest amplitude. 



2.4 Albedo spectrum 

The quantity that we observe is the ratio fp/ f* of the 
spectral-line strengths in the planet's light to those in the 
direct stellar spectrum as a function of wavelength: 



e{a, X) = ^ p(x)gia, X)^ = 6o(A)g(a, A). 



(9) 



This depends on the fraction {Rp/2af of the star's light 
intercepted by the planet, and the fraction of this light that 
is reflected towards the observer. 

The monochromatic stellar flux illuminating the planet, 
orbiting at distance a, is 



t(A) = 



The geometric albedo p{X) is defined at a = 

_ -Fi-oflcctcd(0, A) 



Fin 



t(A) 



(10) 



(11) 



Here -Fi.ofloctod(0, A) = 7r/rcflcctod(0, A) is the disc-averaged 
flux of the starlight reflected from the planet directly back 
toward the star, also measured at the planet's surface. 

As the planet orbits the star, the observed flux varies 
with the star-planet-observer "phase angle" a as the product 
of the geometric albedo p(A) and a "phase function" g{a, A) 
which is normalised to g{0, A) = 1. 

The flux received from the planet at Earth, is therefore 



Ft 

fp{a, A) = p{X)g{a, A)Fincident(A)-^. 



(12) 



Since the stellar flux received at Earth, a distance D 
away, is 



eq. ^ follows. 



(13) 
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Figure 5. A Monte Carlo simulation shows how the gravitational 
acceleration g at the planet's surface is expected to vary with the 
projected orbital velocity amplitude Kp. The system parameters 
are defined as in Section |^. For comparison, gjup = 26.5 m s"'^. 



2.5 Phase function 

The phase angle a depends on the orbital inclination i and 
the orbital phase (p (measured from transit or inferior con- 
junction): 



cos a — ~ sin i. cos 2710. 



(14) 



The observations measure e(a, A) over some range of 
orbital phases <j> and hence phase angles a. However, the 
signal-to-noise ratio and orbital phase coverage of the obser- 
vations is not yet adequate to define the shape of the phase 
function. Accordingly, current practice is to adopt a specific 
phase function in order to express the results in terms of 
the planet/star fiux ratio that would be seen at phase angle 
zero: 



eo(A)=p(A) 



(15) 



Since a is tightly constrained through Kepler's third 
law, given the period P and the star mass M-i, the measure- 
ments of eo(A) measure the product RpyJp^X). 

In this paper we adopt an empirical phase function sim- 
ilar to that of Venus and Jupiter. The expressions given in 
Appendix Dl (eqs. D4 and and Eq. D6 for the phase func- 
tion of a Lambert sphere and Venus respecively) yield the 
phase-dependent flux correction factors plotted in Fig. 0. 



2.6 Planet radius and surface gravity 

Because the surface gravity appears to play an important 
role in determining the height of the cloud deck and hence 
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the albedo, we plot the relationship between g and Kp de- 
rived from the Monte Carlo simulations, in Fig. ^ The 
planet's mass can be determined once the orbital inclina- 
tion is known. Our a priori knowledge of the radius, how- 
ever, comes only from structural models for close-orbiting 
giant planets. 

In these models, the planet radius evolves with time and 
depends on the planet mass. For our purposes, theoretical 
mass-radius-age relations define a range of plausible radii at 
each possible value of the planet mass. The range of po ssible 
planet radii was c omputed specifically for v And b by |GuiL 
lot et 



(1997) 



allowing for uncertainties in the orbital 
inclination and hence the planet's mass. Their radiative- 
convective gas giant models predict a radius between 1.2 
and 1.7 -R,iup for Mp = 0.6Mjup, to between 1.15 and 1.3 
i?jup for Mp = l.lSMjup. 

When this mass-radius relation and its associated un- 
certainty are incorporated in the Monte Carlo model, we 
find that as the inclination decreases, the planet mass in- 
creases and the radius decreases. The lowest values of g are 
therefore found when the orbit is, as is most probable, nearly 
edge-on to the line of sight (Fig. ^). The median predicted 
value, g = 11.6 ms~^, is close to the ~ 10 m s~^ limit below 



which the Cl ass V albedo models of Sudarsky, Burrows & 
Pinto ( L999) ) are applicable. 



If Class V models apply, we can use them to estimate 
the planet's brightness. Viewed fully illuminated, the re- 
flected light from a planet with (Rp/a)^ ~ 10 "'^ and Class 
V albe do p ~ 0.42 (cf . the Class V models of sudarsky. Bur- 
rows & Pinto (1999)) is expected to be 24000 times fainter 



than the direct stellar spectrum in the visual band around 
550 nm. At a phase angle of 60° , however, the reflected spec- 
trum will be a factor two fainter than this (FigQ). 



3 OBSERVATIONS 

The strong dependence of the flux ratio on phase angle indi- 
cates that the best chance of a detection occurs shortly be- 
fore and after superior conjunction. We observed v And with 
the Utrecht Echelle Spectrograph on the 4.2-m William Her- 
schel Telescope at the Roque de los Muchachos Observatory 
on La Palma, on the nights of 2000 Oct 10, Nov 6 and Nov 7. 
These nights were chosen to cover orbital phase ranges when 
the planet is on the far side of the star and well-illuminated, 
and its spectral signature is Doppler-shifted well clear of the 
wings of the stellar profile. The third night was partially af- 
fected by drifting cirrus cloud. Two further nights were also 
lost to cloud. 

The detector was an array of two EEV CCDs, each with 
2048 X 4096 13.5-^im pixels. The CCD was centred at 459.6 
nm in order 124 of the 31 g mm~^ echelle grating, giving 
complete wavelength coverage from 381.3 nm to 675.8 nm 
with minimal vignetting. The average pixel spacing was close 
to 1.6 km s~^, and the full width at half maximum intensity 
of the thorium-argon arc calibration spectra was 3.5 pixels, 
giving an effective resolving power R = 53000. 

Table |^ lists the journal of observations for the three 
nights that contributed to the analysis presented in this pa- 
per. The stellar spectra were exposed for between 300 and 
500 seconds, depending upon seeing, in order to expose the 
CCD to a peak count of 40000 ADU per pixel in the brightest 



parts of the image. A 450-s exposure yielded about 1.2 x 10^ 
electrons per pixel step in wavelength in the brightest or- 
ders in typical (1 arcsec) seeing after extraction. We achieve 
this with the help of an autoguider procedure, which im- 
proves efficiency in good seeing by trailing the stellar image 
up and down the slit by ±2 arcsec during the exposure to 
accumulate the maximum S:N per frame attainable without 
risk of saturation. Note that the 450-s exposure time com- 
pares favourably with the 53-s readout time for the dual 
EEV CCD in terms of observing efficiency - the fraction of 
the time spent collecting photons is above 90 percent. Fol- 
lowing extraction, the S:N in the contiimum of the brightest 
orders is typically 1000 per pixel. 



4 SPECTRUM EXTRACTION 

One-dimensional spectra were extracted from the CCD 
frames using an automated pipeline reduction system built 
around the Starlink ECHOMOP and FIGARO packages. 
Nightly flat-field frames were summed from 50 to 100 frames 
taken at the start and end of each night, using an algo- 
rithm that identifled and rejected cosmic rays and other 
non-repeatable defects by comparing successive frames. The 
nightly flat fields were then added to make a master flat-fleld 
for the entire year's observations. 

The initial tracing of the echelle orders on the CCD 
frames was performed manually on the spectrum of v And 
itself, using exposures taken for this purpose without dither- 
ing the star up and down the slit. The automated extrac- 
tion procedure then subtracted the bias from each frame, 
cropped the frame, determined the form and location of the 
stellar profile on each image relative to the trace, subtracted 
a linear fit to the scattered-light background across the spa- 
tial profile, and performed an optimal (pr o file and inve rse 
variance- weig hted) extraction ( [Horne 198^ [Marsh 1989| ) of 
the orders across the full spatial extent of the object-plus- 
sky region. Flat-field balance factors were applied in the 
process. In all, 62 orders were extracted from each expo- 
sure. The blue CCD recorded orders 148 in the blue to 125 
in the red, giving full spectral coverage from 380.7 to 461.0 
nm with considerable wavelength overlap between adjacent 
orders. Orders 122 to 85 were recorded on the red CCD, cov- 
ering the range 461.9 to 677.3 nm with good overlap. Orders 
123 and 124 fell on the gap between the CCDs, but the loss 
in wavelength coverage was minimal. 



5 STARLIGHT-SUBTRACTED 
VELOCITY-PHASE MAPS 

The maximum expected fiux of starlight scattered from 
V And b is, as we have seen, of order one part in 24000 
of the flux received directly from v And itself. In order to 
detect the planet signal, we flrst subtract the direct stellar 
component from the observed spectrum, leaving the planet 
signal embedded in the residual noise pattern. 

To achieve this, we make a model of the direct starlight 
by aligning and summing all the spectra of the target from 
all nights of the run. The contribution of the planet to this 
summed spectrum is blurred by the planet's orbital motion, 
so that to first order the planet's spectrum is eliminated 
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Table 2. Journal of observations. The UTC mid-times and orbital phases are shown for the first and last groups of four spectra secured 
on each night of observation. The number of groups is given in the final column. 

UTC start Phase UTC End Phase Number 

of groups 

2001 Oct 09 22:27:51 0.296 2001 Oct 10 06:20:26 0.366 18 
2001 Nov 06 19:43:04 0.335 2001 Nov 07 04:27:15 0.413 16 
2001 Nov 07 19:25:48 0.549 2001 Nov 08 04:29:30 0.630 15 



from the composite "template" spectrum. We then scale and 
distort this template spectrum to give the closest possible 
fit to each individual spectrum, using the spline-modulated 
Taylor-series method described in Appendix 

When the aligned template is subtracted from each 
spectrum, we find that the residual spectrum contains a 
spatially fixed but temporally varying pattern of ripples, su- 
perimposed on random noise. The ripples are attributable 
to a combination of time-dependent changes in the detec- 
tor's sensitivity, thermal flexure in the spectrograph causing 
faint ghost reflections to shift slightly on the detector, and 
changes in the strength and velocity of telluric-line absorp- 
tion features in some orders. We map the form of the ripple 
pattern using principal-component analysis, as described in 
Appendix p. After subtracting this map, we are left with 
a planet signal consisting of faint Doppler-shifted copies of 
thousands of stellar absorption lines, deeply buried in noise. 

The positions and identifications of most of these thou- 
sands of lines are well-documented. We use the Vienna 



Phased dynamic spectrum 



atomic-fine database (VALD; Kupka et al. 1999) to compile 
a list of the wavelengths and relative central depths of the 
3450 strongest lines falling in the observed wavelength range, 
for a model atmosphere appropriate to the spectral type and 
surface gravity of the star. We then use least-squares decon- 
volution (LSD; see Appendix ^) to compute a composite 
profile which, when convolved with the line pattern, yields 
an optimal fit to the residual noise spectrum. This procedure 
is similar to cross-correlation in terms of the gain in S:N 
from the weighted summation of thousands of line profiles, 
but has the additional advantage of eliminating sidelobes 
caused by crosstalk with neighbouring lines. 

The resulting devonvolved residual profiles are pre- 
sented in greyscale form as a two-dimensional "velocity- 
phase" map in Fig. ^ This representation of the data shows 
a strong pattern of distortions in the residual stellar profiles 
within ±20 km s^^ of the stellar velocity. These undula- 
tions in the deconvolved profiles appear to be caused by 
high-order mismatches in the spectrum subtraction proce- 
dure. Their amplitude is typically a few parts in 10'' of the 
average continuum level. They vary too rapidly during the 
night to be attributable to, e.g. stellar surface features caus- 
ing time-dependent distortion of the stellar rotation profile. 
A more likely explanation is that the spatial profile produced 
by the telescope dithering procedure was not exactly repeat- 
able from one frame to the next. Given that the slit image 
is slightly tilted with respect to the columns of the detector, 
this could give rise to small changes in the detailed shape of 
the line profile from one exposure to the next. Fortunately 
these ripples only affect a range of velocities at which the 
planet signature would in any case be indistinguishable from 
that of the star. We deliberately avoided observing between 
phases 0.45 and 0.55 for this reason. 
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Figure 6. Velocity-phase map of deconvolved residual profiles de- 
rived from WHT spectra secured on 2000 Oct 09, Nov 06 and Nov 
07. The line weights in the deconvolution were defined assuming 
a grey albedo spectrum. The greyscale runs from black at — 10~* 
times the mean stellar continuum level, to white at -|-10~*. The 
velocity scale is in the reference frame of the star. 



Outside the residual line profile, we would expect to see 
a planet signature as a dark streak following a sinusoidal 
path that crosses the profile from positive to negative ve- 
locity at phase 0.5. No obvious planet signature is, however, 
visible in Fig. ^ 

We verified that a faint planetary signal is preserved 
through the analysis in the presence of realistic noise lev- 
els, by adding a simulated planetary signal to the observed 
spectra and performing the same sequence of operations de- 
scribed above. The simulation procedure simply consisted of 
shifting and scaling the spectrum of v And according to the 
phase function, co-multiplying it by an appropriate geomet- 
ric albedo spectrum, and adding it to the observed data. 

To ensure a strong signal we assumed the planet to have 
a radius 2.0 times that of Jupiter, giving (Rp/a)^ — 2.57 x 
10"*. Initially we used a wavelength-independent geometric 
albedo with p = 0.42, which was chosen to match the contin- 
uum albedo of a Class V model unaffected by alkali-metal 
absorption. When viewed at zero phase angle, the planet- 
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Figure 7. As for Fig. |6| with a simulated planet signature added 
at an orbital inclination of 80°. The model planet has a grey 
albedo spectrum with p = 0.42 and a radius twice that of Jupiter, 
and its signature appears as a dark streak crossing from positive 
to negative velocity at phase 0.5. 
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Figure 8. As for Fig. but showing the matched filter function 
used to measure (Rp/a)-^ at Kp = 137 km s"-"^. The grey scale 
runs from -0.05 (black) to O.Ol(white). Scaling this function to fit 
Fig. ^ (and others like it) yields the strength of the planet signal 
for a given value of Kp. 



to-star flux ratio should thus be eo ~ fp/ f* ~ 1-08 x 10^*. 
We assumed an orbital inclination of 80° , giving an orbital 
velocity amplitude Kp = 137 km and a rotational broad- 
ening of the reflected starlight, Ve,refi = 8 ± 1 km s~^. We 
are therefore justified in using the spectrum of v And itself, 
without any modification of the rotational broadening, to 
approximate the reflected-light spectrum. 

The injected planet signal is clearly visible as a dark 
streak in Fig. |^ Given that no similar signal is easily seen 
in Fig 6, we can conclude that the planet in v And is fainter 
than this one. 



6 EXTRACTING THE PLANET SIGNATURE 

The form of the expected planet signature in the velocity- 
phase diagram can be represented accurately as a travel- 
ling Gaussian of characteristic width Avp whose velocity 
varies sinusoidally around the orbit and whose strength is 
modulated according to the phase function g(a, A) (see Ap- 
pendix Dl). The value of Avp is adjusted to match the ex- 
pected rotational broadening of the reflected starlight. The 



amplitude Kp of the velocity variation and the form of the 
phase function both depend on the orbital inclination i. 

Fig. ^ shows the form of this model signal at an orbital 
inclination i = 80°. The weakening of the simulated plan- 
etary signature near quadrature (phases 0.25 and 0.75) is 
mostly the effect of the phase function. In some cases the 
signal may be further attenuated near quadrature by the 
way in which the templates are computed: since the planet 



signature is nearly stationary in this part of the orbit, some 
of the signal will be removed along with the stellar profile if 
many observations are made in this part of the orbit. The 
planet is detectable on a given night because its velocity 
changes while those of the stellar lines do not. Thus obser- 
vations at quadrature are less helpful. In the present dataset, 
few observations were made near quadrature so this problem 
does not arise. 

The travelling gaussian in this image has nonetheless 
been modified to give an even closer match to the observed 
planet signal, by subtracting its own flux-weighted average. 
This procedure mimics the attenuation of the planet signal 
that occurs when the template spectrum is subtracted from 
the individual spectra. The main effect of the template sub- 
traction is to produce the bright vertical zones seen centred 
at w ~ -70 km s"^ and v ~ -PlOO km s"^ in Fig. 



7 PROBABILITY MAPS FOR Kp AND Rp/A 

We use a sequence of such models for different values of Kp 
as matched filters to measure the strengths and velocity am- 
plitudes of possible faint planet signals of the expected form 
in the velocity-phase maps. At each of a sequence of trial 
values of Kp we construct a model Hij{Kp) like that shown 
in Fig. ^ and scale it to give an optimal fit to the residual 
velocity-phase map Dij using the methods presented in Ap- 
pendix O. For an assumed albedo spectrum p(A) and orbital 
velocity amplitude Kp, the fit of the matched filter to the 
data measures {Rp/af' . 
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Figure 9. The upper panel shows the optimal scaling factor 
(Rp/a)'^ plotted against orbital velocity amplitude Kp, assum- 
ing a grey albedo spectrum with p = 0.42. The lower panel shows 
the associated reduction in measured relative to the fit ob- 
tained in the absence of any planet signal, i.e. for (Rp/a)^ = 0. 
Note that only positive values of {R,p/a)^ are physically plausi- 
ble, so the middle peak in the Ax^ plot is unambiguously a noise 
feature. 
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Figure 10. Relative probability map of model parameters Kp 
and log(€o/p) = log(i?p/a)2, derived from the WHT/UES obser- 
vations of V And, assuming a grey albedo spectrum with p = 0.42. 
The greyscale denotes the probability relative to the best-fit 
model, increasing from for white to 1 for black. From bottom 
to top, the contours show the 1, 2, 3 and 4(7 upper limits limits 
on the signal strength derived from the bootstrap trials The up- 
permost contour, for example, represents the value of log(i?p/a)^ 
that was only exceeded in 3 of 3000 bootstrap trials at each Kp. 
It gives a robust empirical estimate of the upper limit on the 
planet radius allowed by the data for the grey albedo model with 
p = 0.42 assumed here. 



The badness-of-fit statistic gives a relative measure 
of the probabihty that any signal detected could have arisen 
by chance, and is quantified by 



X 



E 



{D,,-{Rp/afH,,{Kp)f 



(16) 



Here afj is the variance associated with Dij, computed from 
the photon statistics of the original image and propagated 
through the deconvolution (see Appendix |c|) . 

In Fig. ^ we plot the optimum scale factor and the as- 
sociated improvement in measured relative to the no- 
planet model. For a noise-dominated signal, the estimated 
value of {Rp/af' will be negative about as often as it is 
positive, and Fig. ^ bears this out. Three candidate peaks 
are seen in the lower panel, of which only the first and the 
third correspond to positive (and therefore physically plau- 
sible) reflected-light fluxes. The third (positive) peak gives 
the greatest improvement in and is therefore the most 
probable, but only by a narrow margin. 

The relative probabilities of the fits to the data for dif- 
ferent values of the free parameters Rp/a and Kp, given by 



P{Kp, [Rp/af) = exp(-x72), 



(17) 



are shown in Fis. |lO| To make this figure, a planet signal pat- 
tern as in Fig. M was computed for many different Kp, and 
each one was scaled by (Rp/a)^ to fit the residuals in Fig ^ 
The greyscale is defined such that white represents the prob- 
ability of a model fit where no planet signal is present, while 
black represents the most probable solution in the map. The 
curves give 1, 2, 3 and 4a detection thresholds derived from 
the bootstrap procedure discussed below. 

To set an upper limit on the strength of the planet sig- 
nal, or to assess the likelihood that a candidate detection is 
spurious, we need to compute the probability of obtaining 
such an improvement in by chance alone. In principle 
this could be done using the distribution for 2 degrees 
of freedom. In practice, however, the distribution of pixel 
values in the deconvolved difference proflles has extended 
non-gaussian tails that demand a more cautious approach. 

Rather than relying solely on formal variances derived 
from photon statistics, we use a "bootstrap" procedure to 
construct empirical distributions for confidence testing, us- 
ing the data themselves. The bootstrap procedure, detailed 
in Appendix E, interchanges at random the order of the 
nights, and rearranges at random the order of spectra in 
each night, thereby scrambling planet signals while retain- 
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Table 3. Upper limits on planet dimensions for various albedo 
models. The upper limits are quoted for an assumed Kp ~ 135 
km s~^, near the peak of the prior probability distribution for 
Kp. Note that the results for the grey albedo model are quoted 
for unit geometric albedo. To obtain the radii for grey models 
with other values of p, the radii in Column 4 should be divided 
by y/P- 
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model 


probability 
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Figure 11. Relative probability map of model parameters Kp 
and log(eo/p) = \og{Rp/a)^ for a simulated planet signature with 
grey albedo p = 0.42, Kp = 137 km s-\ and = 2i?j„p. The 
greyscale and contours are defined as in Fig. hOl The synthetic 
planet signature is detected well above the Act upper limit on the 
strength of noise features in the absence of a planet signal. 



ing the same pattern of systematic errors, phase sampling, 
and statistical noise as in the actual data. 

The 68.4%, 95.4%,99.0% and 99.9% percentage points 
of the resulting bootstrap distribution of (Rp/a)^ at each Kp 
are shown as contours in Fig. ^ and Fig. ^ From bottom 
to top, the contours give the 1, 2, 3, and 4a bootstrap upper 
limits on the strength of the planet signal. They represent 
the signal strengths at which spurious detections occur with 
32, 5, 1, and 0.1% false alarm probability respectively, for 
each fixed value of Kr,. 



7.1 Grey albedo model 

At the most probable values around Kp ~ 135 km s~^, the 
grey albedo model yields a 0.1% bootstrap upper limit on 
the planet/star fiux ratio to < 5.84 x 10~^. Adopting a grey 
albedo model with unit geometric albedo, we find that the 
0.1%, 1.0% and 4.6% upper limits on eo at this velocity 
correspond to upper limits on the planet radius as listed in 
Table |. If the orbital inclination is lower, the planet radius 
is less strongly constrained. 

Two possible planet signals of comparable likelihood are 
seen. Their properties are listed in Table ^. The stronger, at 
Kp ~ 132 km s~^ (i ~ 80°), yields an improvement Ax^ ~ 
13.61 over the model fit obtained assuming no planet signal 
is present. If this feature represents a genuine planet signal, 
its velocity amplitude Kp = 132 km s~^ implies a planet 



mass Mp 



0.74M,jup and (for p 



0.42) a planet radius 
54 



km s~^ (i — 22°) and Ax^ = 11.86, and gives a larger planet 
radius. 

We used the bootstrap simulations to determine the 
probability that a spurious detection with Ax^ > 13.61 
could be produced by a chance alignment of noise features 
in the absence of a genuine planet signal. It is important to 
note that the location of the "blob" between the 2a and 3a 
bootstrap contours in Fig. ^ does not imply a false-alarm 
probability of ~ 3%. These contours give the false-alarm 
probability only if the value of Kp is known in advance, 
which is not the case here. The true false-alarm probability 
is greater, being the frequency with which spurious peaks at 
any plausible value of Kp can exceed the Ax^ of the candi- 
date. If we assume that all values of Kp are equally likely 
in the range 44 km s~^ < Kp < 152 km s~^ over which our 
phase coverage allows signals to be detected, the false-alarm 
probability is found to be 26% via the method described in 
Appendix |e[ 

In practice, however, we are more likely to be fooled 
into believing that a candidate detection near the peak of 
the a priori probability distribution for Kp is genuine, than 
would be the case if the candidate appeared at a velocity 
that was physically implausible given our existing knowl- 
edge of the system parameters. We can therefore refine the 
search range in Kp using our knowledge of the a priori prob- 
ability distribution for Kp, via the method presented in Ap- 
pendix ^. The last two columns of Table ^ show clearly 
that, while the unweighted false alarm probabilities for the 
features near Kp = 132 and Kp — 54 km s~^ are compara- 
ble, the latter is almost certainly spurious. The false-alarm 
probability for the Kp = 132 km s^^ feature drops to 9.4% 
when prior knowledge of Kp is accounted for, while that for 
the 54 km s~^ feature approaches unity. 

For comparison, we show in Fig |ll| a probability map 
derived by subtracting Fig. ^ from Fig. |^ to isolate the 
injected planet signal, then performing the matched-filter 
analysis. The injected planet is clearly detected as a com- 
pact, dark feature with its correct amplitude log{Rp/a)^ = 
2.57 X 10~* and orbital velocity amplitude velocity Kp = 137 
km s"'^. This most probable combination of orbital veloc- 
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Table 4. Orbital velocity amplitude, signal strength, planet radius, Ax^ and false-alarm probabilities (FAP) for possible planet signals. 
The FAP is computed for two different priors. The first prior is uniform in Kp from 44 to 152 km s~^. The second weights Kp in 
proportion to the prior probability based on Call H & K emission and v sin i constraints and the absence of eclipses. In the latter case, 
the false-alarm probability for the Kp = 132 km candidate is found to be between 9 and 10 percent. 
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Figure 12. As for Fig. |^, showing the time-series of decon- 
volved profiles derived from the original observations assuming 
the albedo spectrum to be that of a Class V roaster. 



ity and planet radius yields an improvement Ax^ ~ 98.6 
with respect to the value obtained assuming no planet is 
present (i.e. Rp = 0). This is far greater than the greatest 
Ax^ = 55.3 produced at any value of Kp in any of the boot- 
strap trials. The false-alarm probability for this simulated 
signal is substantially less than one part in 3000, and the 
"detection" is secure. 

We note that the calibration of the (Rp/a)^ scale in 
Fig. ^ depends on the value p = 0.42 assumed for the geo- 
metric albedo. In the next sections, we explore the relative 
ability of non-grey albedo models to fit the data. 
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Figure 13. As for Fig.^, showing the optimal scaling factor 
(Rp/a)^ and the associated improvement in plotted against 
orbital velocity amplitude Kp, assuming a Class V albedo spec- 
trum. 



7.2 Class V roaster model 

The "Class V roaster" is the most highly reflective of the 



models computed by Sudarsky, Burrows & Pinto (1999) 



This model is characteristic of planets with T^jf > 1500 
K and/or surface gravities lower than ~ 10 m s~^, and has 
a silicate cloud deck located high enough in the atmosphere 
that the overlying column density of gaseous alkali metals is 
low, allowing a substantial fraction of incoming photons at 
most optical wavelengths to be scattered back into space. 
There remains, however, a substantial absorption feature 
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Figure 14. Relative probability map of model parameters Kp 
and log(eo/p) = log(i?p/a)2, derived from the WHT/UES obser- 
vations of V And, assuming the albedo spectrum to be that of 
a Class V roaster. The greyscale and contours are defined as in 
Fig- 111- 



around the Na I D lines, as shown in Fig. If v And b 
has a mass and radius close to the values at the peak of the 
prior probability distribution, its surface gravity is close to 
the critical limit (Fig. 

We carried out the deconvolution using the same line list 
as for the grey model, but with the line strengths attenuated 
using the Class V albedo spectrum (see Appendix ^). We 
back-projected the resulting time series of deconvolved pro- 
files (Fig. |l^ as described above. We calibrated the signal 
strength as described in Appendix D4, by injecting an arti- 
ficial planet signature consisting of the spectrum of v And, 
attenuated by the Class V albedo spectrum and scaled to 
the signal strength expected for a planet with Ftp — 2Rjup- 
For the observations we used Avp — 8 km s~^ which again 
yielded the best fit, with scale factors and improvements in 
as plotted in Fig. The probability map for the ob- 
served signals is shown in Fig. hj. 

The form of the Class V probability map is similar to 
that encountered for the grey albedo spectrum. The result- 
ing upper limits on the planet radius are listed in Table 
The local probability maxima near Kp = 52 and 132 km s 
are also present with this albedo model. As in the grey case 



slightly wider margin. Their signal strengths and false-alarm 
probabilities are given in Table ^. 



7.3 Isolated Class IV model 



The Cl ass IV roaster models of ^udarsky, Burrows fc Pintc 
(1999) have a more deeply-buried cloud deck than the Class 
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Figure 15. As for Fig. |6| showing the time-series of decon- 
volved profiles derived from the original observations assuming 
the albedo spectrum to be that of an "isolated" Class IV roaster. 



V models. The resonance lines of Na I and KI are strongly 
saturated, with broad wings due to collisions with H2 ex- 
tending over much of the optical spectrum (Fig. |l]). 

We used the procedures described above to deconvolve 
and back-project the data assuming an "isolated" Class IV 
spectrum. Although this model does not take full account 
of the effects of irradiation of the atmospheric temperature- 
pressure structure, it is a useful compromise between the 
Class V models and the very low albedos found with ir- 
radiated Class IV models. The resulting time-series of de- 
convolved spectra (Fig. |l^) is noisier than the Class V and 
grey-albedo versions, because lines redward of 500 nm con- 
tribute little to the deconvolution. Consequently, the derived 
upper limits on the planet radius (Table H) are higher than 
for the Class V albedo model. 

The plots of {Rp/aY and Ax^ versus Kp show the 
same mix of positive and negative signal amplitudes en- 
countered with the grey and class V spectral models. The 
back-projected probability map (Fig|l7) nonetheless shows 
the Kp = 132 km probability maximum to be present 
and again stronger than the 49 km/sec feature. The fit to 
the data is, however, poorer than in the grey and Class 

V cases, giving an improvement of only Ax^ = 10.8 over 
the no-planet hypothesis and a planetary radius Rp — 
1.68±0.257?j„p (Table §). This is 25% larger than the radius 
estimated with the Class V albedo model. The bootstrap 
test gives a false-alarm probability of 10%. 

This relatively large radius and poor fit appear to arise 
from a mismatch between the wavelength dependence of the 
putative planet signal and the Class IV model. We tested 
this notion by deconvolving and back-projecting a synthetic 
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Figure 16. As for Fig.^, showing tlie optimal scaling factor 
(Rp/a)^ and the associated improvement in plotted against or- 
bital velocity amplitude Kp, assuming an isolated Class IV albedo 
spectrum. 



Class V model using the Class IV line weights and flux cali- 
bration. The radius of the planet was over-estimated by 27% 
in this experiment. This occurs because the least-squares 
deconvolution is attempting to fit a set of lines at all wave- 
lengths, with a deconvolution mask in which the lines at red 
wavelengths are too weak to fit the data properly. The best 
compromise is achieved (in the least-squares sense) by over- 
estimating the depth of the deconvolved profile in order to 
boost the strengths of the weakened lines. 

We conclude that the faint signals that contribute to 
the Kp = 132 km feature are distributed in wavelength 
in a way that resembles more closely a grey or a Class V 
spectrum than a strongly-absorbed Class IV spectrum in 
which only blue light is present. 



7.4 Plausibility of candidate signals 

Comparison of Figs.^, ^ and with the joint prior prob- 
ability distribution in Fig. ^ shows that the feature near 
Kp = 50 km s~^ yields a value of {Rp/af that is consid- 



erably greater than the radius predicted by 



Guillot et al 



(1997), The feature at Kp = 132 ± 3 km s , however, has 



eo = 4.6 X 10 .A plausible geometric albedo p — 0.42 then 
gives {Rp/af = (1.09 ± 0.30) x 10"", which places it very 
close to the peak of the joint prior probability distribution. 

We probed the rotational broadening of the candi- 
date features by performing a series of backprojections us- 
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Figure 17. Relative probability map of model parameters Kp 
and log(eo/p) = log(-Rp/a)^, derived from the WHT/UES obser- 
vations of V And, assuming the albedo spectrum to be that of 
an "isolated" Class IV roaster. The greyscale and contours are 
defined as in Fig. 



ing values of the gaussian width parameter in the range 
6.4 < Avp < 10 km s~^, corresponding to < Ve,refi < 11-3 
km s~^. The feature at 132 km s~^ yields the greatest im- 
provement in when Avp = 8 ± 1 km s~^, or Ue.re/! ~ 
7.0 ± 2 km s~^. The 52 km s~^ feature, however, grows 
less significant as Ve,refi is decreased to the (retrograde) 5 
km s~^ or so expected at the corresponding orbital inclina- 
tion. The rotational broadening of the 132 km s~^ feature 
is therefore consistent with the value predicted in Fig. ^. 

We explored the region of parameter space around both 
features, by carrying out back-projections on a grid of orbital 
periods and epochs of zero phase around the values given in 
Table ^ taking the epoch of transit to be near the middle 
of the data train to ensure that the period and epoch were 
uncorrelated. The local minimum in "^^-s found to be 
centred at Kp = 133 ± 3 km s"\ P = 4.621 ± 0.005 d, and 
To = 2451853.770 ± 0.025. The radial-velocity ephemeris 
gives P = 4.61707 ± 0.00003 d and To = 2451853.777 ± 
0.012, well within the la error bars derived from the back- 
projection. The 52 km s~^ feature, on the other hand, is 
part of a more extended structure whose peak is located 
at JsTp = 67 ± 3 km s'^, P ^ 4.651 ± 0.006 d, and To = 
2451853.81 ± 0.02. This bears little relation to the radial- 
velocity solution, suggesting a spurious origin. 

We conclude that the candidate feature at Kp = 132 
km s~^ corresponds to a local minimum of x^ with respect 
to the parameters of interest, whose rotational broadening, 
phasing and velocity amplitude are consistent with those 
expected of a reflected-light signature from a planet whose 
orbital inclination is near the most probable value. We can- 
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not easily dismiss this feature as being merely an extension 
of a larger, spurious noise-induced structure. 



of PPARC. We thank the referee for raising several useful 
points of clarification. 



8 CONCLUSIONS 

The observations of v And presented here rule out radii for 
the innermost planet greater than Rp > with 
0.1% false-alarm probability if a grey albedo spectrum is as- 
sumed. We derive upper limits also with 0.1% false-alarm 
probability on the planet's radius of 1.53 Rjup assuming the 
albed o spectrum of a Class V roaster (lo w gravity, high cloud 
deck, Sudarsky, Burrows & Pinto 1999), or 2.23 Rjup for an 



isolated Class IV model with saturated Na I D absorption 
from 550 nm to the red limit of our spectra. We cannot, how- 
ever, rule out the possibility that the pl anet has an albedo 
spectrum as b right as a Class V roaster ( Sudarsky, Burrows 
Pinto 199t) if its radius is comparable to that of HD 
209458 b. The evidence for a reflected-light signature in the 
observations reported here, at a projected orbital velocity 
amplitude Kp — 132 ± 3 km s~^, is marginal but encourag- 
ing. If the orbital inclination is as close to edge-on as we infer 
from previously-measured system parameters, the mass and 
radius of the planet yield a surface gravity low enough for 
a Class V atmosphere to be plausible. If future observations 
confirm this signal, it gives a radius Rp = 1.34 ± Q.llRjup 
if a Class V spectrum is assumed. This is comparable to the 
radi us of HD 209458 b inferre d from HST transit photome- 
try (Charbonneau et al. 2000). 

We have explored the spectral properties of the candi- 
date detection via the relative probabilities of the fits to the 
data. The Class V roaster model gives a slightly better fit to 
the data than a grey albedo model. Both yield planet radii 
and orbital velocity amplitudes close to the peak of the prior 
probablility density function. The isolated Class IV albedo 
model gives a significantly poorer fit to the data. 

While the possible detection discussed here gives believ- 
able physical properties for the innermost planet, it remains 
too marginal for us to claim a bona fide detection of refiected 
starlight. Our bootstrap analysis suggests a 7% to 10% like- 
lihood that the feature could be a spurious noise feature. 
There is clearly a strong case for a deeper reflected-light 
search to be made in the v And system. 
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APPENDIX A: REMOVAL OF DIRECT 
STARLIGHT 

In this appendix we describe the procedure we use to sub- 
tract the stellar spectrum without introducing substantial 
noise and also without subtracting the planet signal. The 
intrinsic spectrum of the star was modelled as the co-aligned 
sum of all the spectra of v And secured during the run. 

We began by removing cosmic-ray hits from the indi- 
vidual frames. Within each night's data, each frame was di- 
vided by its predecessor and the result was median smoothed 
with a box covering 21 pixels in the dispersion direction and 
5 adjacent orders. This smoothed frame was used to scale 
the predecessor frame, which was then subtracted from the 
frame under consideration. The difference frame was divided 
by the square root of the frame under consideration, to yield 
deviates in units of the local Poissonian noise amplitude. 
All pixels whose deviates were more than 7 sigma in excess 
of their counterparts in the scaled predecessor frame were 
assumed to be cosmic ray hits. Their values were replaced 
with values taken from the corresponding pixels in the scaled 
predecessor frame. The result was re-scaled and added to 
the scaled predecessor frame to restore the original frame 
with all major cosmic-ray hits removed. The cleaned frames 
were then shifted to co-align the stellar absorption lines to 
the nearest integer pixel. Finally the cleaned fram es w ere 
summed to make the template frame (step 1 in Fig. Al). 



The individual cleaned frames were co-added in groups 
of four for the next stage in the procedure, which involved 
scaling, shifting and blurring (or sharpening) the template 
frame prior to subtraction from each group. The first and 
second derivatives of the template frame were computed 
from the template values tj in each order as a function of 
pixel number Xj in the dispersion direction: 



and 



t 



J + 1 



(Al) 



.(A2) 



The third and fourth derivatives were computed by applying 
the same operations to t" , for use in correcting frame-to- 
frame changes in the higher moments of the PSF and seeing 
profile. 

The vector / of all observed spectrum values within 
each extracted echelle order was then modelled by scaling 
the template t along each order, the scale factors varying as 
a function of position x approximated by a 34-knot least- 
squares spline. The use of such a high-order spline was de- 
manded by vignetting near the edges of each order, which 
produced an abrupt and time-variable change in the slope 
of the stellar continuum towards the end of each order. The 
derivative frames were scaled in a similar fashion, using 6- 
knot splines to define the scale factors. An additional 6-knot 
spline was added to the fit to provide a smooth background 
correction for any inconsistencies in the scattered-light sub- 
traction during the extraction process. The scaled, aligned, 
distorted template spectrum thus had the form 



9j 



J- , /3 -l' I J-" , Z Jll I Jill 



(A3) 



where aj is the value at Xj of the 34-knot spline and Pj, 
7j, 5j, Ej and rjj are the values of the corresponding 6-knot 



splines used to scale the derivatives and correct the back- 
ground. The knot values for the various splines were deter- 
mined by the method of least squares using singular-value 
decomposition to ensure that spurious fiuctuations were sup- 
pressed in those parts of the spectrum where few lines were 
present. 

When the planet velocity is sufficient to shift the ab- 
sorption lines in the refiected starlight well away from those 
in the direct starlight, then the residual spectrum r = f — g 
is expected to retain the reflected starlight signal, Doppler 
shifted and deeply buried in Poissonian noise, once step 2 in 
Fig. Al is complete. 



APPENDIX B: RESIDUAL FIXED-PATTERN 
NOISE 

In this appendix we describe the principal-component anal- 
ysis method used to correct the spatially flxed but time- 
dependent ripples found in the residual spectra following 
subtraction of the template. 

Typically the grouped spectra contained 4 x lO" or more 
photons per pixel over most of the recorded spectrum. The 
scatter of the pixel values in the residual spectrum was com- 
pared with expectations from photon statistics, by dividing 
each grouped spectrum through by the square root of its 
computed variance. The distribution of the pixel deviates 
calculated over the 62 orders used for the deconvolution 
(see Appendix ^ below) was found typically to contain a 
few hundred extreme outliers, mostly caused by systematic 
errors in the polynomial fits in half a dozen regions affected 
by defects on the CCD. Despite clipping at ±4 times the 
expected local root-mean-square (rms) photon noise ampli- 
tude, we found that the distribution of the remaining val- 
ues was Caussian with an rms error between 1.4 and 1.6 
times the expected value. We found a significant correlation 
between pixel values in successive difference frames, indi- 
cating that some fixed-pattern noise sources had not been 
eliminated fully by the template alignment and subtraction. 
The residuals were correlated over a few hours in time. The 
residuals in pairs of frames taken several hours apart or on 
different nights were also found to be correlated, though in 
some instances the slope of the correlation was reduced or 
even reversed. The fixed-pattern noise is visible as ripples in 
the example frame shown following step 2 in Fig. Al. 

We mapped these spatially fixed but time- varying pat- 
terns in the difference spectra using principal-component 
analysis (PCA). We start with a set of M spectra, each 
of length N pixels. The ith spectrum thus has elements 
Si = {sii, Si2, ■ ■ ■ , SiN} and associated variances a^t = 



, o-j]v}- The spatial correlation matrix of the 



time variations in the individual spectral bins of this set 
of spectra is a real symmetric matrix of dimensions N x N, 
whose elements are given by: 



- Sfc 1 



i=l 



Note that we subtract the inverse variance weighted mean 
spectrum s before computing the correlation matrix. The 
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Figure Al. Schematic illustration of the main steps in the construction and subtraction of the template spectrum, and the subsequent 
fixed-pattern noise removal. 



jth spectral bin of s is given by: 

EM 



The principal components are those eigenvectors of the cor- 
relation matrix accounting for the largest fraction of the 
variance, and i n our implementat ion are computed using Ja- 
cobi's method (Press et al. 1992). To keep computing time 
down, we processed the spectra in segments of length 100 
pixels. 

We found that most of the unwanted additional variance 
was contained in the first two principal components, indicat- 
ing that two independent sources of fixed-pattern noise were 
present. The first principal component contained a spatially 
fixed but temporally-varying pattern of ripples afi'ecting all 
orders, which is probably attributable to a slow drift in the 
sensitivity pattern of the fiat field and/or thermal fiexure in 
the spectograph. The second principal component consisted 
mainly of imperfectly-subtracted telluric absorption features 
in some orders. We computed the contributions of these two 
fixed-pattern noise sources to each exposure and subtracted 
them from the difference frames. This effectively removed 
the correlation between successive frames, and reduced the 
RMS scatter in the pixel values to within ±10% of the value 
expected from photon statistics alone (step 3 in Fig. Al). 
The planet signature should not be affected by this proce- 
dure any more than it is affected by the template subtraction 
because, unlike the fixed-pattern noise, it is smeared out by 
orbital motion. 



APPENDIX C: LEAST-SQUARES 
DECONVOLUTION 

We extracted the refiected-light signature from the differ- 
ence frames using the method of weighted least-squar es de- 
convolution (LSD) described by Donati et al. (1997), This 



method, shown schematically in Fig. CI ) entails taking a list 
of lines with relative strengths appropriate to the type of star 
concerned, and computing via the method of least squares a 
"mean" line profile which, when convolved with the line pat- 
tern, gives an optimal match to the observed spectrum. The 
deconvolved profile thus incorporates an average broadening 
function that is representative of all the lines recorded in 
the spectrum. Applied to the residual spectrum from which 
the direct stellar signal has been removed, LSD is an effec- 
tive way of measuring the average line profile of the faint 
refiected-light signature of the planet, because the reflected 
light should should have the same pattern of absorption-line 
positions and strengths as the stellar spectrum. 

In terms of signal improvement, LSD is analogous to 
aligning and averaging (with appropriate weighting factors 
for line strength and local continuum signal strength on 
the recorded frame) the profiles of all the individual pho- 
tospheric absorption lines recorded on each echellogram and 
included in the line list. As each individual spectral line 
appears in at least two adjacent echelle orders, we have 
7700 images of the 3450 spectral lines listed in the observed 
wavelength range. If all lines were of equal strength and the 
continuum signal were constant over the whole frame, the 
signal-to-noise deconvolved profile would be proportional to 
the square root of the number of line images used, i.e. nearly 
70 times greater than the signal-to-noise ratio of a single 
line in the original spectrum. In practice, the recorded con- 
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Figure CI. Schematic illustration of the main steps in the least-squares deconvolution of a time-series of spectra. 



tinuum is not uniform and the lines have a wide variety of 
depths. Even so, the signal-to-noise ratio of the deconvolved 
profile is found to be some 30 times greater than that of 
a single line in the best-exposed parts of the original spec- 
trum. The least-squares fitting procedure has the additional 
advantage that neighbouring, blended lines are treated si- 
multaneously, thereby eliminating the sidelobes that would 
be produced by a simpler shift-and-add procedure. 

We divide the observed residual spectrum by the con- 
tinuum fit_to_yifi_onginal_s2e^ 



residua spectrum r expressed in units of the continuum level 



of the original spectrum. We modify the variances associated 
with the elements of r accordingly. We treat r as the con- 
volution of a "mean" line profile z{v) with a set of weighted 
delta functions at the wavelengths of a comprehensive list 
of spectral lines. The profile z is defined on a linear velocity 
scale, with a velocity increment Ad = 3 km s~^ per bin, 
which is close to the average velocity increment per pixel in 
the extracted spectra. The elements Ajt of the convolution 
matrix A are computed by summing, over all spectral lines 
£, the fractional contribution of the element of the decon- 
volved profile at velocity Vk to the data pixel at wavelength 
Xj when the centre of the deconvolved profile is shifted to 
the wavelength Af of each line in turn. Hence 



(CI) 



The triangular function A has the form A(a;) = 1-1-2; for 
— 1 < a; < 0, A(a;) = 1 ~ x for Q < x < 1, and is zero every- 
where else. The line weights wi, incorporated in the convolu- 
tion matrix A, are proportional to the central depths of the 
lines as computed from a Kurucz model atmosphere for the 
appropriate spectral type. After some experimentation we 



found that a line list computed for a G2 spectral type with 
solar abundances gave the best results. The use of a cooler 
template gives a better fit to the line depths than an earlier 
spectral type, presumably because of the star's above-solar 
metallicity. 

Because the form of the geometric albedo spectrum of 
V And b is unknown, we investigated the effects of both 
grey (i.e. wavelength-independent) and non-grey geometric 
albedo spectra. The non-grey models we em ploved were the 



Class V and "isolate d" Class IV models of judarsky. Bur, 
Pinto (1999)1 (Fig. |). The theoretical albedo spectra 



provided an additional weighting factor in the deconvolution 
linelist, allowing us to place limits on the planet's radius 
for any assumed albedo model, and to assess the relative 
goodness-of-fit to the data for different atmospheric models. 

The deconvolved profile z{v) has the form of a line pro- 
file normalized to unit continuum intensity, but from which 
the continuum has been subtracted. Determination of z{y) 
via least squares entails minimizing the magnitude of the 
misfit vector \r — A ■ z\. 



= {r - A- zf -V ■ {r 



A-z) 



(C2) 



weighted so as to make due allowance for the observational 
errors aj associated with the individual spectral bins. Here 
the inverse variances cr^ associated with the A'^ elements rj 
of the residual spectrum r axe incorporated via the diagonal 
matrix 



y = Diag[l/af, 



, l/cTiv]. 



(C3) 



The least-squares solution for z is found by solving the ma- 
trix equation 



A' - V ■ A-z 



V ■ r. 



(C4) 
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Since the square matrix -V .A is symmetric and positive- 
definite, the least-squares problem can be solve d using effi- 



cien t me thods such as Cholesky decomposition ( Press et al 
199C] 

The deconvolved profile z is expressed in units of the 
weighted mean continuum level. The deconvolution proce- 
dure compensates for local line blends, and so has the ad- 
vantage over cross-correlation methods that outside the re- 
gion occupied by the residual stellar profile, the deconvolved 
spectrum is flat. The formal errors on the M points of the 
deconvolved profile z are obtained in the usual way from the 
diagonal elements of the M x M covariance matrix 



C =\A' - V ■ A] 



(C5) 



planets in our own solar system. Jupiter and Venus appear to 
have phase functions that are more strongly back-scattering 
than a Lambert sphere. Photometric studies of Jupiter at 
large phase a ngles from the Ptoneer and subsequent missions 
have shown (Hovenier & Hage 1989) that the phase function 
for Jupiter is very similar to that of Venus. As a plausible al- 
ternative to the Lambert-sphere formulation, we use a poly- 
nomial approximati on to the em pirically determined phase 
function for Venus (Hilton 1992). The phase-dependent cor- 



rection to the planet's visual magnitude is approximated by: 
Am(Q) = 0.09(Q/100°)+2.39(a/100°)^-0.65(Q/100°)^(D5) 



so that 



10 



-0.4Am(Q) 



(D6) 



APPENDIX D: MATCHED-FILTER ANALYSIS 

To detect the planet we must co-align and add its signal from 
profiles at many different orbit phases. Given an assumed 
orbital inclination, we can compute the sinusoidal path that 
the planet should follow through velocity space, together 
with the attendant changes in signal strength. To search for 
a pattern of faint features displaying the expected behaviour, 
we construct a set of matched filters that can be scaled to 
fit the data for each member of a set of appropriate orbital 



inclinations, as shown schematically in Fig. Dl 



We model the velocity profile of the refiected-light signal 
as a sequence of gaussians with appropriate velocities and 
relative amplitudes: 



G{v,4>\K^) 



exp 



g{<f>,i) X 



1 ( V — Kp sin ( 

AVr, 



(Dl) 



The amplitude Kp of the sinusoidal velocity variation is de- 
termined by the system inclination and stellar mass accord- 
ing to 
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1.3A/f 



1/3 



sin ikm s 



(D2) 



The phase angle a can be determined at any orbital 
phase from 



cos Q = — sin i . cos 2n4). 



(D3) 



Wi, is the integrated area of the deconvolved stellar line 
profile, and Avp is the characteristic width of the planet's 
refiected-light profile. The width Avp = 9.1 km s~^ of the 
gaussian representing the planetary profile was determined 
from a fit to the deconvolved profile of v And. 

Dl Phase function 

For the phase function there are two natural choices. First, 
the phase function of a Lambert sphere is 



g{a) — (sin a + {n ~ a) cos a) /n. 



(D4) 



This assumes that the planetary atmosphere scatters 
isotropically into 27r steradians. 

Second, it may be more realistic to adopt a phase func- 
tion that resembles those for the cloud-covered surfaces of 



D2 Attenuation during starlight subtraction 

The template spectrum that is subtracted in turn from each 
spectrum contains a blurred planet signal, as described in 
Appendix Our matched filter must therefore mimic accu- 
rately the effects of subtracting a template constructed from 
the spectra themselves, so we subtract the inverse variance- 
weighted average of the travelling gaussian. If a^^ is the vari- 
ance of the ith velocity bin in the jth deconvolved profile, 
then the attenuated basis function becomes: 

h,,{Kp) = G{v,,cl,,\Kp) - ' ; ^" \ (D7) 

To allow for any systematic errors in the continuum 
level following template subtraction and deconvolution, we 
subtract the inverse variance-weighted mean value of each 
of the deconvolved profiles at each orbital phase. If dij is the 
original data value at velocity i and phase j, then 



dij 



1/4 



(D8) 



gives the resulting orthogonalised data pixel value. The 
matched filter is orthogonalised in the same way: 



Hi 



hi 



Y.ih^i/^\ 
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(D9) 



D3 Scaling the matched filter 



The attenuated basis function H{v,tp, Kp) is normalised 
such that when it is scaled to give an optimal fit to the 
entire time-series of Dij values, the scaling factor is eg, the 
planet/star flux ratio at a = as deflned in eq. |l^. 

The optimal scale factor eo and its formal error are de- 
termined from 



€0{Kp) = 



D,,H,,{Kp)/al 



^ Hf^{v.,(^,,Kp)/afj 



and 

Var(eo) = J2 



(DIO) 



(Dll) 



We exclude all pixels within 25 km s~^ of the stellar 
line core from the summations, to eliminate spurious effects 
arising from the ripple pattern in the core of the residual 
deconvolved stellar proflle. 
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D4 Calibrating non-grey albedo models. 

The meaning of eo in tlie expressions above is obvious if the 
albedo spectrum is grey, i.e. if eo is independent of wave- 
length. Problems arise, however, when we wish to test how 
well a given non-grey albedo spectrum fits the data. In this 
case, the weighting factors Wi used for the lines in the de- 
convolution mask mask are multiplied by the albedo as de- 
scribed in Appendix ^ above, and eo follows the wavelength 
dependence of the geometric albedo. 

We circumvented this difficult y by replacing eo with 



DIO 



[Rp/a)^ as the scaling factor in eq 

The basis function G is rescaled to become G' — (p) G, 
where (p) is an appropriately weighted wavelength-aver aged 



geometric albedo, then attenuated as before using eq. D7 



For each candidate albedo spectrum we determine the ap- 
propriate value of (p) empirically. We inject an artificial 
planet signal with the required albedo spectrum and known 
Kp and Rp/a into the data, and deconvolve the synthetic 
data with the albedo- weighted line list. We subtract the pro- 
files of the observed spectra, following deconvolution with 
the same line list, to isolate the deconvolved planet signal 
from the noise. We then back-project the isolated planet sig- 
nal, and measure the value of (p) that recovers the correct 
value of Rp/a at the appropriate Kp. 



APPENDIX E: FALSE-ALARM 
PROBABILITIES FOR CANDIDATE SIGNALS 

Candidate reflected-light signals are characterised by posi- 
tive values of (Rp/a)'^ and an improvement in relative 
to the fit obtained when (Rp/a)^ = 0. We determine the 
frequency with which (Rp/a)^ exceeds a given value due to 
noise in the absence of a planet signal, using a bootstrap 
procedure. 

In each of 3000 trials, we randomize the order in which 



the three nights of observations were secured, then we ran- 
domise the order in which the observations were secured 
within each night. The re-ordered observations are then as- 
sociated with the original sequence of dates and times. This 
ensures that any contiguous blocks of spectra containing 
similar systematic errors remain together, but appear at a 
new phase. Any genuine planet signal present in the data 
is, however, completely scrambled in phase. The re-ordered 
data are therefore as capable as the original data of produc- 
ing spurious detections through chance alignments of blocks 
of systematic errors along a single sinusoidal path through 
the data. We record the least-squares estimates of {Rp/af 
and the associated values of as functions of Kp in each 
trial. The 3000 trials implicitly define empirical probability 
distributions of {Rp/aY and that include both the pho- 
ton statistics and the effects of correlated systematic errors 
at each trial value of Kp. 

The percentile points of the distribution of {Rp/af' val- 
ues at each Kp are used to define the la, 2(t, 3(t and 4(t 
upper-limit contours shown in Figs. ^ |l4| and 

To assess the false-alarm probability for a candidate de- 
tection, however, we need to examine the likelihood of the 
fit to the data. In each bootstrap trial, the most likely can- 
didate is the one among those with (Rp/af > that yields 
the greatest improvement Ax^ relative to the no-planet hy- 
pothesis. Such a peak can occur at any value of Kp. We use 
the distribution of the peak Ax'^ value in each trial to deter- 
mine how often we would expect the Ax'^ of the strongest 
spurious noise feature giving (Rp/af > to exceed the ac- 
tual Ax^ of a candidate detection in the observed data. We 
define this probability - summed over all possible values of 
Kp with some weighting according to prior probability esti- 
mates - to be the "false-alarm probability" for a candidate 
signal. 

The unweighted false-alarm probabilities are based on 
the likelihood of obtaining the data D at a given Kp after 
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optimizing (Rp/a)^ , and take no account of whether or not 
that value of Kp is plausible. This conditional likelihood is 
measured relative to the likelihood of getting the same data 
if no planet signal is present: 

PjDlKp) ^ ^^\Kp)-x''('No planet))/2 
P(£)lNo planet) ' ^ ' 

where we have implicitly optimised the value of [Rp/aY 
using the matched-filter fit at each Kp. 

We take our prior knowledge of Kp into account by 
modifying the likelihood to give the joint probability 

P{D, Kp) = P{D\Kp)P{Kp). (E2) 

We use the prior probability distribution for Kp as shown 
projected on to the Kp axis in Fig.^ 

We compute P[D,Kp) for each of the 3000 bootstrap 
trials. To do this, we pick out the greatest value in each trial 
of the quantity 

In P{D, Kp)^\n P{Kp) + ^ (E3) 

where P{Kp) is derived from the simulations described 
in Section 0. The cumulative distribution of maximal 
log P{D , Kp) values derived from the bootstrap trials gives 
the probability that the observed peak value of In P{D, Kp) 
could be spurious. 



